# Library -----------------------------------------------------------------
library(mizer)
## Warning: package 'mizer' was built under R version 4.4.3
library(mizerExperimental)
##
## Attaching package: 'mizerExperimental'
## The following object is masked from 'package:mizer':
##
## validSim
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'tibble' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'stringr' was built under R version 4.4.3
## Warning: package 'forcats' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Setting parameters ------------------------------------------------------
# Set parameters as a single species model
params_initial <- newSingleSpeciesParams()
# Parameters with a reduced resource level
params_rr <- setResource(params_initial,
resource_dynamics = "resource_semichemostat")
# Parameters with initial biomass doubled and initial resources halved
params_double <- params_rr
params_double <- setResource(params_double,
resource_level = 0.1)
initialN(params_double) <- 2*initialN(params_double)
initialNResource(params_double) <- initialNResource(params_double)/2
# Fishing gear ------------------------------------------------------------
# Set up fishing gear using sigmoid_weight() selectivity function
gear_params(params_double) <- data.frame(
gear = "gear",
species = "Target species",
catchability = 0.3,
sel_func = "sigmoid_weight",
sigmoidal_weight = 15,
sigmoidal_sigma = 5)
params_double <- setFishing(params_double, gear_params = gear_params)
gear_params(params_double)
## gear species catchability sel_func
## Target species, gear gear Target species 0.3 sigmoid_weight
## sigmoidal_weight sigmoidal_sigma
## Target species, gear 15 5
# Simulation -------------------------------------------------
# Simulate biomass density when initial biomass is doubled and initial resources
# are reduced
sim_double <- project(params_double, t_max = 15, effort = 1)
animateSpectra(sim_double, total = TRUE, power = 2,
ylim = c(1e-8, NA), wlim = c(1e-3, NA))
# Shows biomass level oscillating --> predator prey relationship
# Extract yield over time dependent on the fishing gear
getYieldGear(sim_double)
## , , sp = Target species
##
## gear
## time gear
## 0 0.0012712274
## 1 0.0004614314
## 2 0.0002662361
## 3 0.0004220103
## 4 0.0008738575
## 5 0.0007496112
## 6 0.0004361817
## 7 0.0003875088
## 8 0.0005553817
## 9 0.0006736737
## 10 0.0005680980
## 11 0.0004720735
## 12 0.0005009610
## 13 0.0005914342
## 14 0.0006062247
## 15 0.0005487976
plotYieldGear(sim_double)
# Figures -----------------------------------------------------------------
# Figure showing flux over increasing biomass density
# Extract weight and count number
N <- finalN(sim_double)["Target species", , drop = TRUE]
w <- w(params_double)
# Calculate individual growth rate
E_growth <- getEGrowth(params_double)["Target species", , drop = TRUE]
gr <- w * E_growth
flux <- gr * N
flux_data <- data.frame(
Weight = w,
Flux = flux
)
# Plot figure
flux_plot <- ggplot(flux_data,
aes(x = Weight, y = Flux)) +
geom_smooth() +
scale_x_continuous(expand = c(0, 0), limits = c(0,105)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
labs(
x = paste0("Weight (g)"),
y = paste0("Flux (g/year)"),
title = "Flux over increasing weight of fish") +
theme_classic()
flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_smooth()`).
ggsave("figures/flux_plot.png",
plot = flux_plot,
device = "png",
width = 8,
height = 6,
units = "in",
dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_smooth()`).
# Plot flux on a log-log axis
flux_plot <- ggplot(flux_data,
aes(x = Weight, y = Flux)) +
geom_smooth() +
scale_x_log10() +
scale_y_log10() +
labs(
x = paste0("Weight (g)"),
y = paste0("Flux (g/year)"),
title = "Flux over increasing weight of fish") +
theme_classic()
flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggsave("figures/flux_plot_with_fishing.png",
plot = flux_plot,
device = "png",
width = 8,
height = 6,
units = "in",
dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# Changing parameters -----------------------------------------------------
# Make local reductions in resource replenishment rate and resource capacity
# Calculate flux
N <- finalN(sim_double)["Target species", , drop = TRUE]
w <- w(params_double)
E_growth <- getEGrowth(params_double)["Target species", , drop = TRUE]
gr <- w * E_growth
flux <- gr * N
# Plot flux before local reductions (without fishing gear)
initial_flux_data <- data.frame(Weight = w,
Flux = flux)
initial_flux_plot <- ggplot(initial_flux_data,
aes(x = Weight, y = Flux)) +
geom_smooth() +
scale_x_log10() +
scale_y_log10() +
labs(
x = paste0("Weight (g)"),
y = paste0("Flux (g/year)"),
title = "Flux over increasing weight of fish") +
theme_classic()
initial_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# Reduce rr
rr <- resource_rate(params_double)
w_full <- w_full(params_double)
w_full[215:240]
## [1] 0.04466836 0.05011872 0.05623413 0.06309573 0.07079458 0.07943282
## [7] 0.08912509 0.10000000 0.11220185 0.12589254 0.14125375 0.15848932
## [13] 0.17782794 0.19952623 0.22387211 0.25118864 0.28183829 0.31622777
## [19] 0.35481339 0.39810717 0.44668359 0.50118723 0.56234133 0.63095734
## [25] 0.70794578 0.79432823
rr[215:240] <- rr[215:240] / 100000000
# Plot resource rate against weight to check reduction
rr_data <- data.frame(
Weight = params_double@w_full,
rr = rr
)
rr_plot <- ggplot(rr_data,
aes(x = Weight, y = rr)) +
geom_smooth() +
scale_x_log10() +
scale_y_log10() +
theme_classic()
rr_plot
## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_smooth()`).
# Reduce resource capacity
rc <- resource_capacity(params_double)
w_full[215:240]
## [1] 0.04466836 0.05011872 0.05623413 0.06309573 0.07079458 0.07943282
## [7] 0.08912509 0.10000000 0.11220185 0.12589254 0.14125375 0.15848932
## [13] 0.17782794 0.19952623 0.22387211 0.25118864 0.28183829 0.31622777
## [19] 0.35481339 0.39810717 0.44668359 0.50118723 0.56234133 0.63095734
## [25] 0.70794578 0.79432823
rc[215:240] <- rc[215:240] / 10000
# Plot resource capacity over weights to see reduction
rc_data <- data.frame(
Weight = params_double@w_full,
rc = rc
)
rc_plot <- ggplot(rc_data,
aes(x = Weight, y = rc)) +
geom_smooth() +
scale_x_log10() +
scale_y_log10() +
theme_classic()
rc_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
params_reduced <- setResource(params_double,
resource_capacity = rc,
resource_rate = rr,
balance = FALSE)
gear_params(params_reduced)
## gear species catchability sel_func
## Target species, gear gear Target species 0.3 sigmoid_weight
## sigmoidal_weight sigmoidal_sigma
## Target species, gear 15 5
# Narrow predation kernel
pred_kernel <- getPredKernel(params_reduced)
pred_kernel_reduced <- pred_kernel[, , 215, drop = FALSE]
ggplot(melt(pred_kernel_reduced)) +
geom_line(aes(x = w_pred, y = value)) +
scale_x_log10()
select(species_params(params_reduced), beta, sigma)
## beta sigma
## Target species 100 1.3
given_species_params(params_reduced)$sigma <- 0.2
given_species_params(params_reduced)$beta <- 100
getPredKernel(params_reduced)[, , 180, drop = FALSE] %>%
melt() %>%
ggplot() +
geom_line(aes(x = w_pred, y = value)) +
scale_x_log10()
# lower maximum intake rate
species_params(params_reduced)$h
## [1] 30
given_species_params(params_reduced)$h <- 30
# Calculate new flux
sim_reduced <- project(params_reduced, t_max = 15, effort = 1)
N_reduced <- finalN(sim_reduced)["Target species", , drop = TRUE]
w <- w(params_reduced)
E_growth_reduced <- getEGrowth(params_reduced)["Target species", , drop = TRUE]
grr <- w * E_growth_reduced
flux_reduced <- grr * N_reduced
reduced_flux_data <- data.frame(Weight = w,
Flux = flux_reduced)
reduced_flux_plot <- ggplot(reduced_flux_data,
aes(x = Weight, y = Flux)) +
geom_smooth() +
scale_x_log10() +
scale_y_log10() +
labs(
x = paste0("Weight (g)"),
y = paste0("Flux (g/year)"),
title = "Flux with locally reduced resource replenishment rate and resource capacity") +
theme_classic()
initial_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
reduced_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# Save plots
ggsave("figures/flux_plot_initial.png",
plot = initial_flux_plot,
device = "png",
width = 8,
height = 6,
units = "in",
dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggsave("figures/flux_plot_reduced.png",
plot = reduced_flux_plot,
device = "png",
width = 8,
height = 6,
units = "in",
dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'